#ifndef MYAMG_RUGESTUBEN_CUH
#define MYAMG_RUGESTUBEN_CUH
#include <Data/CSRMatrixOperator.cuh>
#include <CSRPrecondKernel.cuh>
#include <set>
#include <random>
#include <algorithm>
#include <map>

namespace myamg {
    namespace PointType {
        enum PointTypes {
            UNASSIGNED = 0,
            COARSE = 1,
            FINE = -1,
            ISOLATED = -2,
            S2_FINE = -3
        };
    }

    template<typename Integer, typename Real, typename Allocator>
    class RugeStuben {
    public:
        using matrixType = CSRMatrix<Integer, Real, Allocator>;
        using vectorType = Vector<Real, Allocator>;
        using integerType = Integer;
        using realType = Real;
        using allocatorType = Allocator;

        static void strength(const matrixType &A, Vector<bool, allocatorType> &stg, const double threshold) {
            auto Ap = A.rowPtr();
            auto Ai = A.colIndex();
            auto Av = A.values();
            const auto nRows = A.numRows();
            auto stgData = stg.values();

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                strengthHost(stgData, Ap, Ai, Av, nRows, threshold);
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                dim3 block(512);
                dim3 grid((nRows + block.x - 1) / block.x);
                strengthDevice<<<grid, block>>>(stgData, Ap, Ai, Av, nRows, threshold);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }
        }

        static int HMISCoarsening(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                  Vector<int, allocatorType> &cfmap) {
            const auto numRows = A.numRows();
            const auto numNonzero = A.numNonzeros();
            Integer *Ap;
            Integer *Ai;
            bool *stgData;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                stgData = stg.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[numRows + 1];
                Ai = new Integer[numNonzero];
                stgData = new bool[numNonzero];
                cfmapData = new int[numRows];
                cudaMemcpy(Ap, A.rowPtr(), (numRows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), numNonzero * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(stgData, stg.values(), numNonzero * sizeof(bool), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            auto *weights = new int[numRows];

            auto *stp = new int[numRows + 1];
#pragma omp parallel for
            for (auto i = 0; i <= numRows; i++) {
                stp[i] = 0;
            }

            // 统计每个列的强连接入度
            for (auto i = 0; i < numNonzero; i++) {
                if (stgData[i]) {
                    ++stp[Ai[i]];
                }
            }

            // 累加形成行偏移量
            for (auto i = 0; i < numRows; i++) {
                stp[i + 1] += stp[i];
            }

            const auto nSt = stp[numRows];

            auto *sti = new int[nSt];

            // 填充sti
            for (auto i = static_cast<long>(numRows) - 1; i >= 0; --i) {
                for (auto j = static_cast<long>(Ap[i + 1]) - 1; j >= static_cast<long>(Ap[i]); --j) {
                    if (stgData[j]) {
                        sti[--stp[Ai[j]]] = i;
                    }
                }
            }

            // 计算每个点的权重
            for (auto i = 0; i < numRows; i++) {
                weights[i] = stp[i + 1] - stp[i];
            }

            // 初始化粗点和细点数目
            int nc = 0;
            int nf = 0;

            // 初始化cfmap, 权重为0的点设置为 ISOLATED，其余为 UNASSIGNED
            for (auto i = 0; i < numRows; i++) {
                if (weights[i] == 0) {
                    cfmapData[i] = PointType::ISOLATED;
                    ++nf;
                } else {
                    cfmapData[i] = PointType::UNASSIGNED;
                }
            }

            // 直到所有点被分配为止
            while (nc + nf < numRows) {
                for (auto i = 0; i < numRows; ++i) {
                    auto r = i;
                    // 如果当前点未分配
                    if (cfmapData[i] == PointType::UNASSIGNED) {
                    L0:
                        for (auto j = Ap[r]; j < Ap[r + 1]; ++j) {
                            auto jcol = Ai[j];
                            if (stgData[j] && cfmapData[jcol] == PointType::UNASSIGNED && weights[r] < weights[jcol]) {
                                r = jcol;
                                goto L0;
                            }
                        }

                        for (auto j = stp[r]; j < stp[r + 1]; ++j) {
                            auto jcol = sti[j];
                            if (cfmapData[jcol] == PointType::UNASSIGNED && weights[r] < weights[jcol]) {
                                r = jcol;
                                goto L0;
                            }
                        }

                        cfmapData[r] = PointType::COARSE;
                        ++nc;

                        for (int j = stp[r]; j < stp[r + 1]; ++j) {
                            int jcol = sti[j];
                            if (cfmapData[jcol] == PointType::UNASSIGNED) {
                                cfmapData[jcol] = PointType::FINE;
                                ++nf;

                                for (auto k = Ap[jcol]; k < Ap[jcol + 1]; ++k) {
                                    auto kcol = Ai[k];

                                    if (stgData[k] && cfmapData[kcol] == PointType::UNASSIGNED)
                                        ++weights[kcol];
                                }
                            }
                        }

                        for (auto j = Ap[r]; j < Ap[r + 1]; ++j) {
                            auto jcol = Ai[j];

                            if (stgData[j] && cfmapData[jcol] == PointType::UNASSIGNED) {
                                if (--weights[jcol] == 0) {
                                    cfmapData[jcol] = PointType::FINE;
                                    ++nf;

                                    for (auto k = Ap[jcol]; k < Ap[jcol + 1]; ++k) {
                                        auto kcol = Ai[k];

                                        if (stgData[k] && cfmapData[kcol] == PointType::UNASSIGNED)
                                            ++weights[kcol];
                                    }
                                }
                            }
                        }
                    }
                }
            }

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, numRows * sizeof(int), cudaMemcpyHostToDevice);
                delete[] Ap;
                delete[] Ai;
                delete[] stgData;
                delete[] cfmapData;
            }

            delete[] weights;
            delete[] stp;
            delete[] sti;

            return nc;
        }

        static int HMISLevel2Coarsening(const matrixType &S2, Vector<int, allocatorType> &cfmap) {
            const auto n = S2.numRows();

            Integer *S2p;
            Integer *S2i;
            int *cfmapData;
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                S2p = S2.rowPtr();
                S2i = S2.colIndex();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                S2p = new Integer[n + 1];
                S2i = new Integer[S2.numNonzeros()];
                cfmapData = new int[n];
                cudaMemcpy(S2p, S2.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(S2i, S2.colIndex(), S2.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), n * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
                return -1;
            }

            auto *weight = new Integer[n];

            auto *S2tp = new Integer[n + 1];
            for (auto i = 0; i <= n; ++i) {
                S2tp[i] = 0;
            }

            for (auto j = 0; j < S2.numNonzeros(); ++j) {
                ++S2tp[S2i[j]];
            }

            // 累加形成行偏移量
            for (auto i = 0; i < n; i++) {
                S2tp[i + 1] += S2tp[i];
            }

            // 构建转置的行索引
            auto *S2ti = new Integer[S2tp[n]];
            for (auto i = static_cast<long>(n - 1); i >= 0; --i) {
                for (auto j = static_cast<long>(S2p[i + 1]) - 1; j >= static_cast<long>(S2p[i]); --j) {
                    S2ti[--S2tp[S2i[j]]] = i;
                }
            }

            // 初始化权重weight
            for (auto i = 0; i < n; ++i) {
                weight[i] = S2tp[i + 1] - S2tp[i];
            }

            int nc = 0;
            int nf = 0;
            // 初始化cfmap
            for (auto i = 0; i < n; ++i) {
                cfmapData[i] = PointType::UNASSIGNED;
            }

            while (nc + nf < n) {
                for (auto i = 0; i < n; ++i) {
                    auto r = i;

                    if (cfmapData[r] == PointType::UNASSIGNED) {
                    L0:
                        // 检查行邻居
                        for (auto j = S2p[r]; j < S2p[r + 1]; ++j) {
                            auto jcol = S2i[j];
                            if (cfmapData[jcol] == PointType::UNASSIGNED) {
                                if (weight[r] < weight[jcol]) {
                                    r = jcol;
                                    goto L0;
                                }
                            }
                        }
                        // 检查列邻居
                        for (auto j = S2tp[r]; j < S2tp[r + 1]; ++j) {
                            auto jcol = S2ti[j];
                            if (cfmapData[jcol] == PointType::UNASSIGNED) {
                                if (weight[r] < weight[jcol]) {
                                    r = jcol;
                                    goto L0;
                                }
                            }
                        }

                        // r现在是局部最大权重节点
                        cfmapData[r] = PointType::COARSE;
                        ++nc;

                        // 将r的列邻居标记为细点
                        for (auto j = S2tp[r]; j < S2tp[r + 1]; ++j) {
                            auto jcol = S2ti[j];
                            if (cfmapData[jcol] == PointType::UNASSIGNED) {
                                cfmapData[jcol] = PointType::FINE;
                                ++nf;
                                // 增加细点列邻居的权重
                                for (auto k = S2p[jcol]; k < S2p[jcol + 1]; ++k) {
                                    auto kcol = S2i[k];
                                    if (cfmapData[kcol] == PointType::UNASSIGNED) {
                                        ++weight[kcol];
                                    }
                                }
                            }
                        }

                        // 减少r的行邻居的权重
                        for (auto j = S2p[r]; j < S2p[r + 1]; ++j) {
                            auto jcol = S2i[j];
                            if (cfmapData[jcol] == PointType::UNASSIGNED) {
                                if (--weight[jcol] == 0) {
                                    cfmapData[jcol] = PointType::FINE; // 标记为细点
                                    ++nf;
                                    // 增加jcol行邻居的权重
                                    for (auto k = S2p[jcol]; k < S2p[jcol + 1]; ++k) {
                                        auto kcol = S2i[k];
                                        if (cfmapData[kcol] == PointType::UNASSIGNED) {
                                            ++weight[kcol];
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }

            delete[] weight;
            delete[] S2tp;
            delete[] S2ti;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, n * sizeof(int), cudaMemcpyHostToDevice);
                delete[] S2p;
                delete[] S2i;
                delete[] cfmapData;
            }

            return nc;
        }

        static int AggressiveHMISCoarsening(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                            Vector<int, allocatorType> &cfmap) {
            const auto numRows = A.numRows();

            // 第一级粗化：标准HMIS
            int nc = HMISCoarsening(A, stg, cfmap);

            if (nc == 0) {
                return 0;
            }

            matrixType S2;
            Strength2(A, stg, cfmap, S2);

            // 为第二级粗化准备cfmap
            // aggcf 只包含第一级粗点的标记
            Vector<int, allocatorType> aggcf(nc);

            // 第二级粗化：对粗点再次进行HMIS
            int aggnc = HMISLevel2Coarsening(S2, aggcf);

            // 将两级粗化结果合并
            int *cfmapData;
            int *aggcfData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                cfmapData = cfmap.values();
                aggcfData = aggcf.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cfmapData = new int[numRows];
                aggcfData = new int[nc];
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
                cudaMemcpy(aggcfData, aggcf.values(), nc * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
                return -1;
            }

            // 合并两级结果
            for (int i = 0, j = 0; i < numRows; ++i) {
                if (cfmapData[i] >= 0) {
                    // 如果是第一级的粗点
                    cfmapData[i] = aggcfData[j++]; // 用第二级的结果替换
                }
                // 如果是细点或孤立点，保持不变
            }

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, numRows * sizeof(int), cudaMemcpyHostToDevice);
                delete[] cfmapData;
                delete[] aggcfData;
            }

            return aggnc;
        }

        static int PMISCoarsening(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                  Vector<int, allocatorType> &cfmap) {
            const auto numRows = A.numRows();
            const auto numNonzero = A.numNonzeros();

            Integer *Ap;
            Integer *Ai;
            bool *stgData;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                stgData = stg.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[numRows + 1];
                Ai = new Integer[numNonzero];
                stgData = new bool[numNonzero];
                cfmapData = new int[numRows];
                cudaMemcpy(Ap, A.rowPtr(), (numRows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), numNonzero * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(stgData, stg.values(), numNonzero * sizeof(bool), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
                return -1;
            }

            const int lambda = 32768;
            auto *weights = new int[numRows];

            // 初始化基础权重（0-32767）
            std::random_device rd;
            std::mt19937 gen(rd());
            std::uniform_int_distribution<int> dist(0, 32767);

            for (int i = 0; i < numRows; ++i) {
                weights[i] = dist(gen);
            }

            // 计算每个点的强连接入度
            // weights[i] += lambda * (强连接到 i 的边数)
            for (int j = 0; j < numNonzero; ++j) {
                if (stgData[j]) {
                    weights[Ai[j]] += lambda;
                }
            }

            int nc = 0;
            int nf = 0;

            // 初始化 cfmap：权重 < lambda 的点标记为 ISOLATED
            for (int i = 0; i < numRows; ++i) {
                if (weights[i] < lambda) {
                    cfmapData[i] = PointType::ISOLATED;
                    ++nf;
                } else {
                    cfmapData[i] = PointType::UNASSIGNED;
                }
            }

            int LESS = 0;
            int GREATER = 2;

            // 迭代直到所有点被分配
            while (nc + nf < numRows) {
                std::swap(LESS, GREATER);

                for (int i = 0; i < numRows; ++i) {
                    if (cfmapData[i] == GREATER) {
                        for (int j = Ap[i]; j < Ap[i + 1]; ++j) {
                            int jcol = Ai[j];

                            if (stgData[j] && (cfmapData[jcol] == GREATER || cfmapData[jcol] == LESS)) {
                                if (weights[jcol] > weights[i] ||
                                    (weights[jcol] == weights[i] && jcol < i)) {
                                    cfmapData[i] = LESS;
                                    break;
                                } else {
                                    cfmapData[jcol] = LESS;
                                }
                            }
                        }
                    }
                }

                for (int i = 0; i < numRows; ++i) {
                    if (cfmapData[i] == GREATER) {
                        cfmapData[i] = PointType::COARSE;
                        ++nc;
                    }
                }

                for (int i = 0; i < numRows; ++i) {
                    if (cfmapData[i] == LESS) {
                        for (int j = Ap[i]; j < Ap[i + 1]; ++j) {
                            if (stgData[j] && cfmapData[Ai[j]] == PointType::COARSE) {
                                cfmapData[i] = PointType::FINE;
                                ++nf;
                                break;
                            }
                        }
                    }
                }
            }

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, numRows * sizeof(int), cudaMemcpyHostToDevice);
                delete[] Ap;
                delete[] Ai;
                delete[] stgData;
                delete[] cfmapData;
            }

            delete[] weights;

            return nc;
        }

        static int AggressivePMISCoarsening(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                            Vector<int, allocatorType> &cfmap) {
            const auto numRows = A.numRows();

            // 第一级粗化：标准 PMIS
            int nc = PMISCoarsening(A, stg, cfmap);

            if (nc == 0) {
                return 0;
            }

            // 构造二阶强连接矩阵 S2
            matrixType S2;
            Vector<bool, allocatorType> stg2;
            Strength2(A, stg, cfmap, S2);

            Vector<int, allocatorType> aggcf(nc);

            // 第二级粗化：对粗点再次进行 PMIS
            int aggnc = PMISCoarsening(S2, stg2, aggcf);

            // 将两级粗化结果合并
            int *cfmapData;
            int *aggcfData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                cfmapData = cfmap.values();
                aggcfData = aggcf.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cfmapData = new int[numRows];
                aggcfData = new int[nc];
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
                cudaMemcpy(aggcfData, aggcf.values(), nc * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
                return -1;
            }

            for (int i = 0, j = 0; i < numRows; ++i) {
                if (cfmapData[i] >= 0) {
                    cfmapData[i] = aggcfData[j++];
                }
            }

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, numRows * sizeof(int), cudaMemcpyHostToDevice);
                delete[] cfmapData;
                delete[] aggcfData;
            }

            return aggnc;
        }

        static void Strength2(const matrixType &A, const Vector<bool, allocatorType> &stg,
                              const Vector<int, allocatorType> &cfmap,
                              matrixType &S2) {
            const auto numRows = A.numRows();
            const auto numNonzero = A.numNonzeros();

            Integer *Ap;
            Integer *Ai;
            bool *stgData;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                stgData = stg.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[numRows + 1];
                Ai = new Integer[numNonzero];
                stgData = new bool[numNonzero];
                cfmapData = new int[numRows];
                cudaMemcpy(Ap, A.rowPtr(), (numRows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), numNonzero * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(stgData, stg.values(), numNonzero * sizeof(bool), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
                return;
            }

            // 重新编号粗点，并统计粗点数量
            int nc = 0;
            auto *coarseMap = new int[numRows];
            for (int i = 0; i < numRows; ++i) {
                if (cfmapData[i] >= 0) {
                    coarseMap[i] = nc++;
                } else {
                    coarseMap[i] = -1;
                }
            }

            // 第二步：构造 Sc
            auto *Sc_rowptr = new Integer[numRows + 1];
            Sc_rowptr[0] = 0;

            for (int i = 0; i < numRows; ++i) {
                if (cfmapData[i] < 0) {
                    // i 是细点或孤立点
                    int cnt = 0;
                    for (int j = Ap[i]; j < Ap[i + 1]; ++j) {
                        if (stgData[j] && cfmapData[Ai[j]] >= 0) {
                            ++cnt;
                        }
                    }
                    Sc_rowptr[i + 1] = cnt;
                } else {
                    Sc_rowptr[i + 1] = 0;
                }
            }

            for (int i = 0; i < numRows; ++i) {
                Sc_rowptr[i + 1] += Sc_rowptr[i];
            }

            auto *Sc_colind = new Integer[Sc_rowptr[numRows]];

            for (int i = 0; i < numRows; ++i) {
                if (cfmapData[i] < 0) {
                    // i 是细点或孤立点
                    int k = Sc_rowptr[i];
                    for (int j = Ap[i]; j < Ap[i + 1]; ++j) {
                        if (stgData[j] && cfmapData[Ai[j]] >= 0) {
                            // 存储映射后的粗点索引
                            Sc_colind[k++] = coarseMap[Ai[j]];
                        }
                    }
                }
            }

            auto *S2_rowptr = new Integer[nc + 1];
            S2_rowptr[0] = 0;

            // 使用工作数组来标记访问过的列，避免重复
            auto *marker = new int[nc];
            for (int i = 0; i < nc; ++i) {
                marker[i] = -1;
            }

            // 先统计每行的非零元个数
            for (int i = 0; i < numRows; ++i) {
                if (cfmapData[i] >= 0) {
                    // i 是粗点
                    int ci = coarseMap[i];
                    int cnt = 0;

                    // 遍历 i 的所有强连接邻居 k
                    for (int kk = Ap[i]; kk < Ap[i + 1]; ++kk) {
                        if (!stgData[kk]) continue;

                        int k = Ai[kk];

                        // 遍历 k 的所有粗点邻居 j（通过 Sc）
                        for (int jj = Sc_rowptr[k]; jj < Sc_rowptr[k + 1]; ++jj) {
                            int cj = Sc_colind[jj];

                            // 排除对角元，并避免重复计数
                            if (cj != ci && marker[cj] != ci) {
                                marker[cj] = ci;
                                ++cnt;
                            }
                        }
                    }

                    S2_rowptr[ci + 1] = cnt;
                }
            }

            // 累加形成行偏移
            for (int i = 0; i < nc; ++i) {
                S2_rowptr[i + 1] += S2_rowptr[i];
            }

            // 分配列索引数组
            Integer nnz = S2_rowptr[nc];
            auto *S2_colind = new Integer[nnz];

            // 重置标记数组
            for (int i = 0; i < nc; ++i) {
                marker[i] = -1;
            }

            // 填充列索引
            for (int i = 0; i < numRows; ++i) {
                if (cfmapData[i] >= 0) {
                    // i 是粗点
                    int ci = coarseMap[i];
                    int idx = S2_rowptr[ci];

                    // 遍历 i 的所有强连接邻居 k
                    for (int kk = Ap[i]; kk < Ap[i + 1]; ++kk) {
                        if (!stgData[kk]) continue;

                        int k = Ai[kk];

                        // 遍历 k 的所有粗点邻居 j（通过 Sc）
                        for (int jj = Sc_rowptr[k]; jj < Sc_rowptr[k + 1]; ++jj) {
                            int cj = Sc_colind[jj];

                            // 排除对角元，并避免重复
                            if (cj != ci && marker[cj] != ci) {
                                marker[cj] = ci;
                                S2_colind[idx++] = cj;
                            }
                        }
                    }
                }
            }

            // 清理 Sc
            delete[] Sc_rowptr;
            delete[] Sc_colind;
            delete[] marker;
            delete[] coarseMap;

            // 将结果拷贝到 S2 矩阵
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                S2 = matrixType(nc, nc, nnz);

                std::copy(S2_rowptr, S2_rowptr + nc + 1, S2.rowPtr());
                std::copy(S2_colind, S2_colind + nnz, S2.colIndex());

                delete[] S2_rowptr;
                delete[] S2_colind;
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                S2 = matrixType(nc, nc, nnz);

                cudaMemcpy(S2.rowPtr(), S2_rowptr, (nc + 1) * sizeof(Integer), cudaMemcpyHostToDevice);
                cudaMemcpy(S2.colIndex(), S2_colind, nnz * sizeof(Integer), cudaMemcpyHostToDevice);

                delete[] S2_rowptr;
                delete[] S2_colind;
            }

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Ap;
                delete[] Ai;
                delete[] stgData;
                delete[] cfmapData;
            }
        }

        static void longRangeInterpolation(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                           Vector<int, allocatorType> &cfmap,
                                           matrixType &P, const int minElements, const int maxElements,
                                           const double truncationFactor) {
            auto numRows = A.numRows();

            directInterpolation(A, stg, cfmap, P);

            auto *mask = new bool[numRows];
            Integer *Ap;
            Integer *Pp;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Pp = P.rowPtr();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[numRows + 1];
                Pp = new Integer[numRows + 1];
                cfmapData = new int[numRows];
                cudaMemcpy(Ap, A.rowPtr(), (numRows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pp, P.rowPtr(), (numRows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            for (int i = 0; i < numRows; ++i) {
                auto numP = Pp[i + 1] - Pp[i];
                auto numN = Ap[i + 1] - Ap[i] - 1;
                mask[i] = cfmapData[i] >= 0 || numP >= minElements || numN == numP;
            }

            delete[] cfmapData;
            delete[] Ap;
            delete[] Pp;

            jacobiInterpolation(A, P, mask);
            delete[] mask;
            interpTruncation(P, cfmap, maxElements, truncationFactor);
        }

        static void aggressiveLongRangeInterpolation(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                                     Vector<int, allocatorType> &cfmap,
                                                     matrixType &P, const int minElements, const int maxElements,
                                                     const double truncationFactor) {
            auto n = A.numRows();

            directInterpolation(A, stg, cfmap, P);

            auto *mask = new bool[n];
            Integer *Pp;
            int *cfmapData;
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Pp = P.rowPtr();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Pp = new Integer[n + 1];
                cfmapData = new int[n];
                cudaMemcpy(Pp, P.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), n * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            for (auto i = 0; i < n; ++i) {
                mask[i] = cfmapData[i] >= 0 || (cfmapData[i] == PointType::FINE && (Pp[i + 1] > Pp[i]));
            }

            jacobiInterpolation(A, P, mask);

            for (auto i = 0; i < n; ++i) {
                mask[i] = cfmapData[i] >= 0 || cfmapData[i] == PointType::S2_FINE;
            }

            jacobiInterpolation(A, P, mask);

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Pp = P.rowPtr();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Pp;
                Pp = new Integer[n + 1];
                cudaMemcpy(Pp, P.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            for (auto i = 0; i < n; ++i) {
                mask[i] = cfmapData[i] >= 0 || (Pp[i + 1] - Pp[i]) >= minElements;
            }
            jacobiInterpolation(A, P, mask);

            delete[] mask;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Pp;
                delete[] cfmapData;
            }

            interpTruncation(P, cfmap, maxElements, truncationFactor);
        }

        static void directInterpolation(const matrixType &A, const Vector<bool, allocatorType> &stg,
                                        Vector<int, allocatorType> &cfmap,
                                        matrixType &P) {
            auto numRows = A.numRows();
            Integer *Ap;
            Integer *Ai;
            Real *Av;
            bool *stgData;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
                stgData = stg.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[numRows + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                stgData = new bool[A.numNonzeros()];
                cfmapData = new int[numRows];
                cudaMemcpy(Ap, A.rowPtr(), (numRows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(stgData, stg.values(), A.numNonzeros() * sizeof(bool), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), numRows * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            // 将COARSE点的cfmap值重新映射为编号
            auto nc = 0;
            for (auto i = 0; i < numRows; ++i) {
                if (cfmapData[i] >= 0) {
                    // COARSE点
                    cfmapData[i] = nc++;
                }
            }

            auto *Prowptr = new Integer[numRows + 1];
            Prowptr[0] = 0;

            for (auto i = 0; i < numRows; ++i) {
                if (cfmapData[i] < 0) {
                    // FINE点或者ISOLATED点
                    int cnt = 0;
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        auto jcol = Ai[j];
                        if (stgData[j] && cfmapData[jcol] >= 0) {
                            // 强连接且为COARSE点
                            cnt++;
                        }
                    }
                    Prowptr[i + 1] = cnt;
                } else {
                    // COARSE点
                    Prowptr[i + 1] = 1;
                }
            }

            // 累加形成行偏移量
            for (auto i = 0; i < numRows; ++i) {
                Prowptr[i + 1] += Prowptr[i];
            }

            // 分配P的其余空间
            auto *PcolIndex = new Integer[Prowptr[numRows]];
            auto *Pvalues = new Real[Prowptr[numRows]];

            for (auto i = 0; i < numRows; ++i) {
                if (cfmapData[i] < 0) {
                    // FINE点或者ISOLATED点
                    Real diag = 0;
                    Real sumN = 0; // 非对角元素之和
                    Real sumP = 0; // 强连接的COARSE点权重和

                    // 找对角元和计算sumN
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        if (Ai[j] == i) {
                            diag = Av[j];
                        } else {
                            sumN += Av[j];
                        }
                    }

                    // 填充P矩阵，计算sumP
                    auto k = Prowptr[i];
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        auto jcol = Ai[j];
                        if (stgData[j] && cfmapData[jcol] >= 0) {
                            // 强连接且为COARSE点
                            auto val = Av[j];
                            auto jval = ((val > 0) ^ (diag > 0)) ? val : static_cast<Real>(0);
                            PcolIndex[k] = cfmapData[jcol];
                            Pvalues[k++] = jval;
                            sumP += jval;
                        }
                    }

                    Real omega = 0;
                    if (sumP != 0 && diag != 0) {
                        omega = -sumN / (diag * sumP);
                    }

                    for (auto j = Prowptr[i]; j < Prowptr[i + 1]; ++j) {
                        Pvalues[j] *= omega;
                    }
                } else {
                    // COARSE点
                    PcolIndex[Prowptr[i]] = cfmapData[i];
                    Pvalues[Prowptr[i]] = 1.0;
                }
            }

            // 构造矩阵P
            matrixType tempP(numRows, nc, Prowptr[numRows]);
            tempP.setRowPtr(Prowptr);
            tempP.setColIndex(PcolIndex);
            tempP.setValues(Pvalues);
            P = tempP;
            delete[] Prowptr;
            delete[] PcolIndex;
            delete[] Pvalues;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, numRows * sizeof(int), cudaMemcpyHostToDevice);
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
                delete[] stgData;
                delete[] cfmapData;
            }
        }

        static void jacobiInterpolation(const matrixType &A, matrixType &P, const bool *mask) {
            auto n = A.numRows(); // A的行数
            auto m = P.numCols(); // P的列数

            Integer *Ap;
            Integer *Ai;
            Real *Av;
            Integer *Pp;
            Integer *Pi;
            Real *Pv;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
                Pp = P.rowPtr();
                Pi = P.colIndex();
                Pv = P.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[n + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                Pp = new Integer[n + 1];
                Pi = new Integer[P.numNonzeros()];
                Pv = new Real[P.numNonzeros()];
                cudaMemcpy(Ap, A.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pp, P.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pi, P.colIndex(), P.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pv, P.values(), P.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            // 统计Q矩阵每一行的非零元个数
            auto *Qrowptr = new Integer[n + 1];
            Qrowptr[0] = 0;

            // 使用标记数组避免重复计数
            auto *w = new int[m];
            for (auto i = 0; i < m; ++i) {
                w[i] = -1;
            }

            for (auto i = 0; i < n; ++i) {
                if (mask[i]) {
                    Qrowptr[i + 1] = Pp[i + 1] - Pp[i];
                } else {
                    int cnt = 0;
                    for (auto k = Ap[i]; k < Ap[i + 1]; ++k) {
                        auto kcol = Ai[k];
                        if (kcol == i) {
                            continue;
                        }
                        for (auto j = Pp[kcol]; j < Pp[kcol + 1]; ++j) {
                            auto jcol = Pi[j];
                            if (w[jcol] != i) {
                                w[jcol] = i;
                                cnt++;
                            }
                        }
                    }
                    Qrowptr[i + 1] = cnt;
                }
            }
            // 累加形成行偏移量
            for (auto i = 0; i < n; ++i) {
                Qrowptr[i + 1] += Qrowptr[i];
            }

            // 分配Q的其余空间
            auto *QcolIndex = new Integer[Qrowptr[n]];
            auto *Qvalues = new Real[Qrowptr[n]];

            // 重置标记数组
            for (auto i = 0; i < m; ++i) {
                w[i] = -1;
            }

            for (auto i = 0; i < n; ++i) {
                if (mask[i]) {
                    for (auto j = Pp[i], k = Qrowptr[i]; j < Pp[i + 1]; ++j, ++k) {
                        QcolIndex[k] = Pi[j];
                        Qvalues[k] = Pv[j];
                    }
                } else {
                    // 计算 Q[i] = A[i] * P
                    auto r = Qrowptr[i];
                    auto r0 = r;
                    Real diag = 0;
                    Real sumN = 0; // A的第i行所有非对角元之和
                    Real sumQ = 0; // Q的第i行所有元素之和

                    // 找到A的对角元
                    for (auto k = Ap[i]; k < Ap[i + 1]; ++k) {
                        if (Ai[k] == i) {
                            diag = Av[k];
                            break;
                        }
                    }

                    // 计算 Q[i] = A[i] * P
                    for (auto k = Ap[i]; k < Ap[i + 1]; ++k) {
                        auto kcol = Ai[k];
                        if (kcol == i) {
                            continue;
                        }
                        auto kval = Av[k];
                        auto _kval = ((kval > 0) ^ (diag > 0)) ? kval : static_cast<Real>(0);
                        sumN += kval;

                        for (auto j = Pp[kcol]; j < Pp[kcol + 1]; ++j) {
                            auto jcol = Pi[j];
                            auto jval = Pv[j];
                            if (w[jcol] < r0) {
                                w[jcol] = r;
                                QcolIndex[r] = jcol;
                                Qvalues[r++] = _kval * jval;
                            } else {
                                Qvalues[w[jcol]] += _kval * jval;
                            }
                        }
                    }

                    // 计算sumQ
                    for (auto j = Qrowptr[i]; j < Qrowptr[i + 1]; ++j) {
                        sumQ += Qvalues[j];
                    }

                    // 计算并应用缩放因子omega
                    Real omega = 0;
                    if (sumQ != 0 && diag != 0) {
                        omega = -sumN / (diag * sumQ);
                    }
                    for (auto j = Qrowptr[i]; j < Qrowptr[i + 1]; ++j) {
                        Qvalues[j] *= omega;
                    }
                }
            }

            delete[] w;

            // 构造矩阵Q
            matrixType Q(n, m, Qrowptr[n]);
            Q.setRowPtr(Qrowptr);
            Q.setColIndex(QcolIndex);
            Q.setValues(Qvalues);

            delete[] Qrowptr;
            delete[] QcolIndex;
            delete[] Qvalues;

            // 用Q替换P
            P = Q;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
                delete[] Pp;
                delete[] Pi;
                delete[] Pv;
            }
        }

        static void interpTruncation(matrixType &P, Vector<int, allocatorType> &cfmap, const int maxElements,
                                     const double truncationFactor) {
            auto n = P.numRows(); // P的行数
            auto nc = P.numCols(); // P的列数
            Integer *Pp;
            Integer *Pi;
            Real *Pv;
            int *cfmapData;
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Pp = P.rowPtr();
                Pi = P.colIndex();
                Pv = P.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Pp = new Integer[n + 1];
                Pi = new Integer[P.numNonzeros()];
                Pv = new Real[P.numNonzeros()];
                cfmapData = new int[n];
                cudaMemcpy(Pp, P.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pi, P.colIndex(), P.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pv, P.values(), P.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), n * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            auto *Qrowptr = new Integer[n + 1];
            Qrowptr[0] = 0;

            // 处理每一行，进行截断和重新缩放
            for (auto i = 0; i < n; ++i) {
                if (cfmapData[i] < 0) {
                    // FINE点
                    auto m = Pp[i + 1] - Pp[i]; // 当前行非零元个数
                    if (m == 0) {
                        Qrowptr[i + 1] = 0;
                    } else {
                        std::vector<std::pair<Integer, Real> > entries(m);
                        // 复制该行的数据
                        for (auto j = Pp[i], k = 0; j < Pp[i + 1]; ++j, ++k) {
                            entries[k] = {Pi[j], Pv[j]};
                        }
                        std::sort(entries.begin(), entries.end(),
                                  [](const std::pair<Integer, Real> &a, const std::pair<Integer, Real> &b) -> bool {
                                      return std::abs(a.second) > std::abs(b.second);
                                  });

                        auto theta = truncationFactor * std::abs(static_cast<double>(entries[0].second));
                        Real sum = 0; // 保留元素的和
                        Real drop = 0; // 丢弃元素的和
                        int cnt = 0; // 保留元素的计数

                        for (auto j = 0; j < m; ++j) {
                            if (std::abs(entries[j].second) < theta || j == maxElements) {
                                break;
                            }
                            Pi[Pp[i] + cnt] = entries[j].first;
                            Pv[Pp[i] + cnt] = entries[j].second;
                            sum += entries[j].second;
                            cnt++;
                        }

                        for (auto j = cnt; j < m; ++j) {
                            drop += entries[j].second;
                        }

                        // 重新缩放保留的元素
                        Real omega = 0;
                        if (sum != 0) {
                            omega = (sum + drop) / sum;
                        }
                        for (auto j = Pp[i]; j < Pp[i] + cnt; ++j) {
                            Pv[j] *= omega;
                        }
                        Qrowptr[i + 1] = cnt;
                    }
                } else {
                    // COARSE点
                    Qrowptr[i + 1] = 1;
                }
            }

            // 累加形成行偏移量
            for (auto i = 0; i < n; ++i) {
                Qrowptr[i + 1] += Qrowptr[i];
            }

            // 将截断后的数据恢复到新的数组中
            auto *QcolIndex = new Integer[Qrowptr[n]];
            auto *Qvalues = new Real[Qrowptr[n]];

            for (auto i = 0; i < n; ++i) {
                for (auto j = Pp[i], k = Qrowptr[i]; k < Qrowptr[i + 1]; ++j, ++k) {
                    QcolIndex[k] = Pi[j];
                    Qvalues[k] = Pv[j];
                }
            }

            // 构造新的P矩阵
            matrixType Q(n, nc, Qrowptr[n]);
            Q.setRowPtr(Qrowptr);
            Q.setColIndex(QcolIndex);
            Q.setValues(Qvalues);
            P = Q;
            delete[] Qrowptr;
            delete[] QcolIndex;
            delete[] Qvalues;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                cudaMemcpy(cfmap.values(), cfmapData, n * sizeof(int), cudaMemcpyHostToDevice);
                delete[] Pp;
                delete[] Pi;
                delete[] Pv;
                delete[] cfmapData;
            }
        }

        static void smoothedAggregation(const matrixType &A, const Vector<bool, allocatorType> stg,
                                        Vector<int, allocatorType> &cfmap, matrixType &P) {
            auto n = A.numRows();

            Integer *Ap;
            Integer *Ai;
            Real *Av;
            bool *stgData;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
                stgData = stg.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[n + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                stgData = new bool[A.numNonzeros()];
                cfmapData = new int[n];
                cudaMemcpy(Ap, A.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(stgData, stg.values(), A.numNonzeros() * sizeof(bool), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), n * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            // 将COARSE点的cfmap值重新映射为编号
            int nc = 0;
            for (auto i = 0; i < n; ++i) {
                if (cfmapData[i] >= 0) {
                    cfmapData[i] = nc++;
                }
            }

            // 为每个FINE点分配聚合
            auto *aggregates = new int[n];
            for (auto i = 0; i < n; ++i) {
                if (cfmapData[i] < 0) {
                    // FINE点
                    int iagg = -1; // 表示未分配
                    Real amax = -1;
                    // 找到强连接中绝对值最大的COARSE点
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        if (!stgData[j]) {
                            continue;
                        }
                        auto jcol = Ai[j];
                        auto absVal = std::abs(Av[j]);
                        if (cfmapData[jcol] >= 0 && absVal > amax) {
                            iagg = cfmapData[jcol];
                            amax = absVal;
                        }
                    }
                    aggregates[i] = iagg;
                } else {
                    aggregates[i] = cfmapData[i];
                }
            }

            // 处理未分配的FINE点
            auto *mask = new bool[n];
            for (auto i = 0; i < n; ++i) {
                mask[i] = aggregates[i] != -1;
            }

            // 为未分配的点找到权重最大的邻居聚合
            for (auto i = 0; i < n; ++i) {
                if (!mask[i]) {
                    std::map<int, Real> iaggwgt; // 聚合索引到权重的映射
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        auto jcol = Ai[j];
                        auto absVal = std::abs(Av[j]);
                        if (mask[jcol]) {
                            iaggwgt[aggregates[jcol]] += absVal;
                        }
                    }

                    // 选择权重最大的聚合
                    int iagg = -1;
                    Real amax = -1;
                    for (auto iter = iaggwgt.begin(); iter != iaggwgt.end(); ++iter) {
                        if (iter->second > amax) {
                            iagg = iter->first;
                            amax = iter->second;
                        }
                    }
                    aggregates[i] = iagg;
                }
            }

            delete[] mask;

            // 计算对角元的倒数
            auto *recipD = new Real[n];
            for (auto i = 0; i < n; ++i) {
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    if (Ai[j] == i) {
                        recipD[i] = 1.0 / Av[j];
                        break;
                    }
                }
            }

            // 计算阻尼因子omega
            Real omega = 4.0 / (3.3 * largestEigenvalueEstimate(A, recipD, 10));

            // 统计P矩阵每行非零元个数
            auto *Prowptr = new Integer[n + 1];
            Prowptr[0] = 0;

            auto *w = new int[nc];
            for (auto i = 0; i < nc; ++i) {
                w[i] = -1;
            }

            for (auto i = 0; i < n; ++i) {
                int cnt = 0;
                if (aggregates[i] != -1) {
                    w[aggregates[i]] = i;
                    ++cnt;
                }
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    auto jcol = Ai[j];
                    if (jcol == i) {
                        continue;
                    }
                    auto jagg = aggregates[jcol];
                    if (jagg != -1 && w[jagg] != i) {
                        w[jagg] = i;
                        ++cnt;
                    }
                }
                Prowptr[i + 1] = cnt;
            }

            // 累加形成行偏移量
            for (auto i = 0; i < n; ++i) {
                Prowptr[i + 1] += Prowptr[i];
            }

            // 填充P矩阵
            auto PcolIndex = new Integer[Prowptr[n]];
            auto Pvalues = new Real[Prowptr[n]];

            // 重置标记数组
            for (auto i = 0; i < nc; ++i) {
                w[i] = -1;
            }

            for (auto i = 0; i < n; ++i) {
                auto k = Prowptr[i];
                auto k0 = k;

                // 添加自己的聚合
                if (aggregates[i] != -1) {
                    w[aggregates[i]] = k;
                    PcolIndex[k] = aggregates[i];
                    Pvalues[k++] = 1.0 - omega;
                }

                auto multiplier = -omega * recipD[i];

                // 添加邻居的聚合
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    auto jcol = Ai[j];
                    if (jcol == i) {
                        continue;
                    }
                    auto jval = Av[j];
                    auto jagg = aggregates[jcol];

                    if (jagg == -1) {
                        continue;
                    }

                    if (w[jagg] < k0) {
                        w[jagg] = k;
                        PcolIndex[k] = jagg;
                        Pvalues[k++] = multiplier * jval;
                    } else {
                        Pvalues[w[jagg]] += multiplier * jval;
                    }
                }
            }

            delete[] w;
            delete[] aggregates;
            delete[] recipD;

            matrixType temp(n, nc, Prowptr[n]);
            temp.setRowPtr(Prowptr);
            temp.setColIndex(PcolIndex);
            temp.setValues(Pvalues);
            P = temp;
            delete[] Prowptr;
            delete[] PcolIndex;
            delete[] Pvalues;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
                delete[] stgData;
                delete[] cfmapData;
            }
        }

        static void aggressiveSmoothedAggregation(const matrixType &A, const Vector<bool, allocatorType> stg,
                                                  Vector<int, allocatorType> &cfmap, matrixType &P) {
            auto n = A.numRows();

            Integer *Ap;
            Integer *Ai;
            Real *Av;
            bool *stgData;
            int *cfmapData;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
                stgData = stg.values();
                cfmapData = cfmap.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[n + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                stgData = new bool[A.numNonzeros()];
                cfmapData = new int[n];
                cudaMemcpy(Ap, A.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(stgData, stg.values(), A.numNonzeros() * sizeof(bool), cudaMemcpyDeviceToHost);
                cudaMemcpy(cfmapData, cfmap.values(), n * sizeof(int), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            // 将COARSE点的cfmap值重新映射为编号
            int nc = 0;
            for (auto i = 0; i < n; ++i) {
                if (cfmapData[i] >= 0) {
                    cfmapData[i] = nc++;
                }
            }

            // 为每个FINE点分配聚合
            auto *aggregates = new int[n];
            for (auto i = 0; i < n; ++i) {
                if (cfmapData[i] < 0) {
                    // FINE点
                    int iagg = -1; // 表示未分配
                    Real amax = -1;
                    // 找到强连接中绝对值最大的COARSE点
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        if (!stgData[j]) {
                            continue;
                        }
                        auto jcol = Ai[j];
                        auto absVal = std::abs(Av[j]);
                        if (cfmapData[jcol] >= 0 && absVal > amax) {
                            iagg = cfmapData[jcol];
                            amax = absVal;
                        }
                    }
                    aggregates[i] = iagg;
                } else {
                    aggregates[i] = cfmapData[i];
                }
            }

            // 三次聚合传递
            auto *mask = new bool[n];
            for (auto k = 0; k < 3; ++k) {
                // 标记已分配聚合的点
                for (auto i = 0; i < n; ++i) {
                    mask[i] = aggregates[i] != -1;
                }
                // 为未分配的点找到权重最大的邻居聚合
                for (auto i = 0; i < n; ++i) {
                    if (!mask[i]) {
                        std::map<int, Real> iaggwgt; // 聚合索引到权重的映射
                        for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                            auto jcol = Ai[j];
                            auto absVal = std::abs(Av[j]);
                            if (mask[jcol]) {
                                iaggwgt[aggregates[jcol]] += absVal;
                            }
                        }

                        // 选择权重最大的聚合
                        int iagg = -1;
                        Real amax = -1;
                        for (auto iter = iaggwgt.begin(); iter != iaggwgt.end(); ++iter) {
                            if (iter->second > amax) {
                                iagg = iter->first;
                                amax = iter->second;
                            }
                        }
                        aggregates[i] = iagg;
                    }
                }
            }

            delete[] mask;

            auto *recipD = new Real[n];

            for (auto i = 0; i < n; ++i) {
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    if (Ai[j] == i) {
                        recipD[i] = 1.0 / Av[j];
                        break;
                    }
                }
            }

            // 计算阻尼因子omega
            Real omega = 4.0 / (3.3 * largestEigenvalueEstimate(A, recipD, 10));

            // 统计P矩阵每行非零元个数
            auto *Prowptr = new Integer[n + 1];
            Prowptr[0] = 0;

            auto *w = new int[nc];
            for (auto i = 0; i < nc; ++i) {
                w[i] = -1;
            }

            for (auto i = 0; i < n; ++i) {
                int cnt = 0;
                if (aggregates[i] != -1) {
                    w[aggregates[i]] = i;
                    ++cnt;
                }
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    auto jcol = Ai[j];
                    if (jcol == i) {
                        continue;
                    }
                    auto jagg = aggregates[jcol];
                    if (jagg != -1 && w[jagg] != i) {
                        w[jagg] = i;
                        ++cnt;
                    }
                }
                Prowptr[i + 1] = cnt;
            }

            // 累加形成行偏移量
            for (auto i = 0; i < n; ++i) {
                Prowptr[i + 1] += Prowptr[i];
            }

            // 填充P矩阵
            auto PcolIndex = new Integer[Prowptr[n]];
            auto Pvalues = new Real[Prowptr[n]];

            // 重置标记数组
            for (auto i = 0; i < nc; ++i) {
                w[i] = -1;
            }

            for (auto i = 0; i < n; ++i) {
                auto k = Prowptr[i];
                auto k0 = k;

                // 添加自己的聚合
                if (aggregates[i] != -1) {
                    w[aggregates[i]] = k;
                    PcolIndex[k] = aggregates[i];
                    Pvalues[k++] = 1.0 - omega;
                }

                auto multiplier = -omega * recipD[i];

                // 添加邻居的聚合
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    auto jcol = Ai[j];
                    if (jcol == i) {
                        continue;
                    }
                    auto jval = Av[j];
                    auto jagg = aggregates[jcol];

                    if (jagg == -1) {
                        continue;
                    }

                    if (w[jagg] < k0) {
                        w[jagg] = k;
                        PcolIndex[k] = jagg;
                        Pvalues[k++] = multiplier * jval;
                    } else {
                        Pvalues[w[jagg]] += multiplier * jval;
                    }
                }
            }

            delete[] w;
            delete[] aggregates;
            delete[] recipD;

            // 构造P矩阵
            matrixType temp(n, nc, Prowptr[n]);
            temp.setRowPtr(Prowptr);
            temp.setColIndex(PcolIndex);
            temp.setValues(Pvalues);
            P = temp;
            delete[] Prowptr;
            delete[] PcolIndex;
            delete[] Pvalues;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
                delete[] stgData;
                delete[] cfmapData;
            }
        }

        static Real largestEigenvalueEstimate(const matrixType &A, const Real *recipD, const int numIteration) {
            auto n = A.numRows();
            Integer *Ap;
            Integer *Ai;
            Real *Av;
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[n + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                cudaMemcpy(Ap, A.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }
            auto *x = new Real[n];
            auto *y = new Real[n];
            // 初始化x为全1向量
            for (auto i = 0; i < n; ++i) {
                x[i] = 1.0;
            }
            // 幂迭代
            for (auto iter = 0; iter < numIteration; ++iter) {
                // y = D^(-1) * A * x
                for (auto i = 0; i < n; ++i) {
                    Real sum = 0;
                    for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                        sum += Av[j] * x[Ai[j]];
                    }
                    y[i] = recipD[i] * sum;
                }

                // 归一化
                Real norm = 0;
                for (auto i = 0; i < n; ++i) {
                    norm += y[i] * y[i];
                }
                norm = std::sqrt(norm);
                for (auto i = 0; i < n; ++i) {
                    x[i] = y[i] / norm;
                }
            }

            // 计算Rayleigh商，lambda = x^T * D^(-1) * A * x
            Real lambda = 0;
            for (auto i = 0; i < n; ++i) {
                Real sum = 0;
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    sum += Av[j] * x[Ai[j]];
                }
                lambda += x[i] * recipD[i] * sum;
            }

            delete[] x;
            delete[] y;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
            }

            return lambda;
        }

        static void RAP(
            const matrixType &R,
            const matrixType &A,
            const matrixType &P,
            matrixType &C) {
            auto R_rows = R.numRows();
            auto A_rows = A.numRows();
            auto A_cols = A.numCols();
            auto P_cols = P.numCols();

            Integer *Rp, *Ri, *Ap, *Ai, *Pp, *Pi;
            Real *Rv, *Av, *Pv;

            // 根据内存类型获取数据指针
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Rp = R.rowPtr();
                Ri = R.colIndex();
                Rv = R.values();
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
                Pp = P.rowPtr();
                Pi = P.colIndex();
                Pv = P.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Rp = new Integer[R_rows + 1];
                Ri = new Integer[R.numNonzeros()];
                Rv = new Real[R.numNonzeros()];
                Ap = new Integer[A_rows + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                Pp = new Integer[A_cols + 1];
                Pi = new Integer[P.numNonzeros()];
                Pv = new Real[P.numNonzeros()];

                cudaMemcpy(Rp, R.rowPtr(), (R_rows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ri, R.colIndex(), R.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Rv, R.values(), R.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ap, A.rowPtr(), (A_rows + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pp, P.rowPtr(), (A_cols + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pi, P.colIndex(), P.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Pv, P.values(), P.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
                return;
            }

            // 第一步: 计算 AP = A * P
            std::vector<std::map<Integer, Real> > AP_temp(A_rows);

            for (Integer i = 0; i < A_rows; ++i) {
                for (Integer idx_a = Ap[i]; idx_a < Ap[i + 1]; ++idx_a) {
                    Integer k = Ai[idx_a];
                    Real a_val = Av[idx_a];

                    for (Integer idx_p = Pp[k]; idx_p < Pp[k + 1]; ++idx_p) {
                        Integer j = Pi[idx_p];
                        Real p_val = Pv[idx_p];
                        AP_temp[i][j] += a_val * p_val;
                    }
                }
            }

            // 构造 AP 的 CSR 格式
            auto *AProwptr = new Integer[A_rows + 1];
            AProwptr[0] = 0;
            for (Integer i = 0; i < A_rows; ++i) {
                AProwptr[i + 1] = AProwptr[i] + AP_temp[i].size();
            }

            Integer AP_nnz = AProwptr[A_rows];
            auto *APcolIndex = new Integer[AP_nnz];
            auto *APvalues = new Real[AP_nnz];

            Integer idx = 0;
            for (Integer i = 0; i < A_rows; ++i) {
                for (const auto &entry: AP_temp[i]) {
                    APcolIndex[idx] = entry.first;
                    APvalues[idx] = entry.second;
                    idx++;
                }
            }

            // 第二步: 计算 C = R * AP
            std::vector<std::map<Integer, Real> > C_temp(R_rows);

            for (Integer i = 0; i < R_rows; ++i) {
                for (Integer idx_r = Rp[i]; idx_r < Rp[i + 1]; ++idx_r) {
                    Integer k = Ri[idx_r];
                    Real r_val = Rv[idx_r];

                    for (Integer idx_ap = AProwptr[k]; idx_ap < AProwptr[k + 1]; ++idx_ap) {
                        Integer j = APcolIndex[idx_ap];
                        Real ap_val = APvalues[idx_ap];
                        C_temp[i][j] += r_val * ap_val;
                    }
                }
            }

            // 清理 AP 临时数据
            delete[] AProwptr;
            delete[] APcolIndex;
            delete[] APvalues;

            // 第三步: 构造 C 矩阵的 CSR 格式
            auto *Crowptr = new Integer[R_rows + 1];
            Crowptr[0] = 0;
            for (Integer i = 0; i < R_rows; ++i) {
                Crowptr[i + 1] = Crowptr[i] + C_temp[i].size();
            }

            Integer C_nnz = Crowptr[R_rows];
            auto *CcolIndex = new Integer[C_nnz];
            auto *Cvalues = new Real[C_nnz];

            idx = 0;
            for (Integer i = 0; i < R_rows; ++i) {
                for (const auto &entry: C_temp[i]) {
                    CcolIndex[idx] = entry.first;
                    Cvalues[idx] = entry.second;
                    idx++;
                }
            }

            // 构造新的 C 矩阵
            matrixType temp(R_rows, P_cols, C_nnz);
            temp.setRowPtr(Crowptr);
            temp.setColIndex(CcolIndex);
            temp.setValues(Cvalues);
            C = temp;

            // 清理内存
            delete[] Crowptr;
            delete[] CcolIndex;
            delete[] Cvalues;

            // 如果是 Cuda 内存，需要释放临时分配的 Host 内存
            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Rp;
                delete[] Ri;
                delete[] Rv;
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
                delete[] Pp;
                delete[] Pi;
                delete[] Pv;
            }
        }

        static void sparsification(const matrixType& A, matrixType& B, const double threshold) {
            auto n = A.numRows();
            Integer *Ap;
            Integer *Ai;
            Real *Av;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Host) {
                Ap = A.rowPtr();
                Ai = A.colIndex();
                Av = A.values();
            } else if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                Ap = new Integer[n + 1];
                Ai = new Integer[A.numNonzeros()];
                Av = new Real[A.numNonzeros()];
                cudaMemcpy(Ap, A.rowPtr(), (n + 1) * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Ai, A.colIndex(), A.numNonzeros() * sizeof(Integer), cudaMemcpyDeviceToHost);
                cudaMemcpy(Av, A.values(), A.numNonzeros() * sizeof(Real), cudaMemcpyDeviceToHost);
            } else {
                std::cerr << "Unknown allocator type" << std::endl;
            }

            // 每行和每列的非对角元素最大绝对值
            auto* rmax = new Real[n];
            auto* cmax = new Real[n];

            // 初始化
            for (auto i = 0; i < n; ++i) {
                rmax[i] = 0;
                cmax[i] = 0;
            }

            // 计算行最大值和列最大值
            for (auto i = 0; i < n; ++i) {
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    auto jcol = Ai[j];
                    if (jcol == i) {
                        continue;
                    }

                    auto absJval = std::abs(Av[j]);
                    if (absJval > rmax[i]) {
                        rmax[i] = absJval;
                    }
                    if (absJval > cmax[jcol]) {
                        cmax[jcol] = absJval;
                    }
                }
            }

            // 标记哪些元素需要保留
            auto* keep = new bool[Ap[n]];
            auto* Browptr = new Integer[n + 1];
            Browptr[0] = 0;

            // 判断哪些元素需要保留，统计每行的非零元素数
            for (auto i = 0; i < n; ++i) {
                int cnt = 1; // 对角元素始终保留
                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    auto jcol = Ai[j];
                    if (jcol == i) {
                        keep[j] = false; // 对角元素稍后处理
                    } else {
                        auto absJval = std::abs(Av[j]);
                        keep[j] = (absJval > rmax[i] * static_cast<Real>(threshold)) || (absJval > cmax[jcol] * static_cast<Real>(threshold));
                        if (keep[j]) {
                            cnt++;
                        }
                    }
                }
                Browptr[i + 1] = cnt;
            }

            // 累加形成行偏移量
            for (auto i = 0; i < n; ++i) {
                Browptr[i + 1] += Browptr[i];
            }

            // 分配新矩阵的空间
            auto* BcolIndex = new Integer[Browptr[n]];
            auto* Bvalues = new Real[Browptr[n]];

            // 填充新矩阵，并计算对角线补偿
            for (auto i = 0; i < n; ++i) {
                Real diag = 0;
                auto k = Browptr[i] + 1; // 跳过对角线元素的位置

                for (auto j = Ap[i]; j < Ap[i + 1]; ++j) {
                    if (keep[j]) {
                        BcolIndex[k] = Ai[j];
                        Bvalues[k] = Av[j];
                        ++k;
                    } else {
                        diag += Av[j];
                    }
                }
                // 设置对角线元素
                BcolIndex[Browptr[i]] = i;
                Bvalues[Browptr[i]] = diag;
            }

            delete[] rmax;
            delete[] cmax;
            delete[] keep;

            matrixType temp(n, n, Browptr[n]);
            temp.setRowPtr(Browptr);
            temp.setColIndex(BcolIndex);
            temp.setValues(Bvalues);
            delete[] Browptr;
            delete[] BcolIndex;
            delete[] Bvalues;
            B = temp;

            if constexpr (AllocatorType<allocatorType>::value == MemoryType::Cuda) {
                delete[] Ap;
                delete[] Ai;
                delete[] Av;
            }
        }
    };
}

#endif //MYAMG_RUGESTUBEN_CUH
